On the Occurrence of Finite-Time-Singularities in Epidemic Models of Rupture, 

Earthquakes and Starquakes 
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We present a new kind of critical stochastic finite-time-singularity, relying on the interplay be- 
tween long-memory and extreme fluctuations. We illustrate it on the well-established epidemic-type 
aftershock (ETAS) model for aftershocks, based solely on the most solidly documented stylized 
facts of seismicity (clustering in space and in time and power law Gutenberg-Richter distribution 
of earthquake energies). This theory accounts for the main observations (power law acceleration 
and discrete scale invariant structure) of critical rupture of heterogeneous materials, of the largest 
sequence of starquakes ever attributed to a neutron star as well as of earthquake sequences. 
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A large portion of the current work on rupture and 
earthquake prediction is based on the search for precur- 
sors to large events in the seismicity itself. Observations 
of the acceleration of seismic moment leading up to large 
events and "stress shadows" following them have been 
interpreted as evidence that seismic cycles represent the 
approach to and retreat from a critical state of a fault 
network pj. This "critical state" concept is fundamen- 
tally different from the long-time view of the crust as 
evolving spontaneously in a statistically stationary crit- 
ical state, called self-organized criticality (SOC) |2|. In 
the SOC view, all events belong to the same global popu- 
lation and participate in shaping the self-organized criti- 
cal state. Large earthquakes are inherently unpredictable 
because a big earthquake is simply a small earthquake 
that did not stop. By contrast, in the critical point view, 
a great earthquake plays a special role and signals the end 
of a cycle on its fault network. The dynamical organiza- 
tion is not statistically stationary but evolves as the great 
earthquake becomes more probable. Predictability might 
then become possible by monitoring the approach of the 
fault network towards the critical state. This hypothesis 
first proposed in Q] is the theoretical induction of a se- 
ries of observations of accelerated seismicity |3|, Li which 
has been later strengthened by several other observations 
[M El 13 • Theoretical support has also come from sim- 
ple computer models of critical rupture |9j and experi- 
ments of material rup ture , cellular automata, with 
[l ll| and without Il2| long-range interaction, and from 
granular simulators Models of regional seismicity 

with more faithful fault geometry have been developed 
that also show accelerating seismicity before large model 
events 0HOH. 

There are at least five different mechanisms that are 
known to lead to critical accelerated seismicity of the 
form 



N(t) oc i/(t c -ty 



(1) 



ending at the critical time t c , where N(t) is the seismic- 
ity rate (or acoustic emission rate for material rupture). 
Such finite-time-singularities are quite common and have 
been found in many well-established models of natural 
systems, either at special points in space such as in the 
Euler equations of inviscid fluids, in vortex collapse of 
systems of point vortices, in the equations of General Rel- 
ativity coupled to a mass field leading to the formation 
of black holes, in models of micro-organisms aggregating 
to form fruiting bodies, or in the more prosaic rotating 
coin (Euler 's disk). They all involve some kind of positive 
feedback, which in the rupture context can be the follow- 
ing (see ^3 f° r a review) : sub-critical crack growth , 
geometrical feedback in creep rupture feedback of 
damage on the elastic coefficients with strain dependent 
damage rate feedback in a percolation model of re- 
gional seismicity 

.17], 

feedback in a stress-shadow model 
for regional seismicity |15l Il7| . 

While these mechanisms are plausible, their relevance 
to the earth crust remains unproven. Here, we present 
a novel mechanism leading to a new kind of critical 
stochastic finite-time-singularity in the seismicity rate, 
using the well-established epidemic-type aftershock se- 
quence (ETAS) model for aftershocks, introduced by 
[2fJ,|21|, based solely on the most solidly documented styl- 
ized facts of seismicity mentioned above. The adjective 
"stochastic" emphasizes the fact that the critical time t c 
is determined in large part by the specific sets of inno- 
vations of the random process. We show that, in a finite 
domain of its parameter space, the rate of seimic activity 
in the ETAS model diverges in finite time according to 
(fT|) . The underlying mechanism relies on large deviations 
occuring in an explosive branching process. One of the 
advantage of this discovery is to be able to account for 
the observations of accelerated seismicity and acoustic 
emission in material failure, without invoking any new 
ingredient other than those already well-established em- 
pirically. We apply this insight to quantify the longest 
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available starquake sequence of a neutron star soft 7-ray 
repeaters. 

We shall use the example of earthquakes but the model 
applies similarly to microcracking in materials. The 
ETAS model is a generalization of the modified Omori 
law, in that it takes into account the secondary after- 
shock sequences triggered by all events. The modified 
Omori's law states that the occurrence rate of the di- 
rect aftershock-daughters from an earthquake decreases 
with the time from the mainshock according to the "bare 
propagator" K/(t + c) p . In the ETAS model, all earth- 
quakes are simultaneously mainshocks, aftershocks and 
possibly foreshocks. Contrary to the usual definition of 
aftershocks, the ETAS model does not impose an after- 
shock to have an energy smaller than the mainshock. 
This way, the same law describes both foreshocks, af- 
tershocks and mainshocks. An observed "aftershock" se- 
quence of a given earthquake (starting the clock) is the 
result of the activity of all events triggering events trig- 
gering themselves other events, and so on, taken together. 
The corresponding seismicity rate (the "dressed propa- 
gator"), which is given by the superposition of the after- 
shock sequences of all events, is the quantity we derive 
here. 

Each earthquake (the "mother") of energy Ei > Eq 
occurring at time gives birth to other events ( "daugh- 
ters" ) of any possible energy, chosen with the Gutenberg- 
Richter density distribution P(E) = fi/(E/E ) 1+ ^ with 
exponent [i ~ 2/3, at a later time between t and t + dt 
at the rate 

Mt-U) =p(Ei) ¥(*-*<) ■ (2) 

p{Ei) — K (Ei/E ) a gives the number of daughters 
born from a mother with energy Ei, with the same ex- 
ponent a for all earthquakes. This term accounts for 
the fact that large mothers have many more daughters 
than small mothers because the larger spatial exten- 
sion of their rupture triggers a larger domain. Eq is a 
lower bound energy below which no daughter is triggered. 
ty(t— ti) = ( t _i + c }i+ti is the normalized waiting time dis- 
tribution (local Omori's law or "bare propagator") giving 
the rate of daughters born a time t — t, after the mother. 

The ETAS model is fundamentally a "branching" 
model |2-l| with no "loops", i.e., each event has a unique 
"mothcr-mainshock" and not several. This "mean-field" 
or random phase approximation allows us to simplify the 
analysis while still keeping the essential physics in a qual- 
itative way. The problem is to calculate the "dressed" or 
"renormalized" propagator (rate of seismic activity) that 
includes the whole cascade of secondary sequences [25| . 
The key parameter is the average number n (or "branch- 
ing ratio" ) of daughter-earthquakes created per mother- 
event, summed over all possible energies, n is equal to 
the integral of 4>Ei(t — ti) over all times after U and over 
all energies Ei > Eq. This integral converges to a finite 



value n < 00 for 8 > (local Omori's law decay faster 
than 1 /t) and for a < \i (not too large a growth of the 
number of daughters as a function of the energy of the 
mother). The resulting average rate N(t) of seismicity is 
the solution of the Master equation |26( 

N{t)= / dr N(t) / dE' P(E') <p E ,(t-r) (3) 

Jo Je 

giving the number N(t)dt of events occurring between t 
and t + dt of any possible energy. We have made explicit 
the upper bound E max (t) equal to the typical maximum 
earthquake energy sampled up to time t. For a < /1, 
this upper bound has no impact on the results and can 
be replaced by +00 |26j. There may be a source term 
S(t) to add to the r.h.s. of ©, corresponding to cither 
a constant background seismicity or to a large triggering 
earthquake. In this last case, the rate N(t) solution of 
© is the "dressed" propagator giving the renormalized 
Omori's law. A rich behavior, which has been fully clas- 
sified by a complete analytical treatment |26|. has been 
found: sub-criticality n < 1 |25| and super-criticality 
n > 1 26], where n depends on the control parameters 
/1, a, 8, K and c. With a single value of the exponent 
1 + 8 of the "bare" propagator ^>(t) ~ l/t 1+e , we obtain 
a continuum of apparent exponents for the global rate 
of aftershocks [26j which may account for the observed 
variability of Omori's exponent p around p — 1 reported 
by many workers. 

Here, we explore the regime a > /i, for which n is in- 
finite. This signals the impact of large earthquake ener- 
gies, suggesting the relevance of the upper bound E max (t) 
in • This case is actually observed in real seismicity by 
|27j ] who obtained a > fj, for some aftershock sequences 
in Greece, and by [3] who found a > /1 for 13 out of 34 
aftershock sequences in Japan. This case a > (i also char- 
acterizes the seismic activity preceding the 1984 M = 6.8 
Nagano Prefecture earthquake |22| . After the mainshock, 
the seismicity returned in the sub-critical regime 8 > 0, 
a < \l and n < 1. 

This case a > fi is similar to that found underlying var- 
ious situations of anomalous transport 28] : in this regime 
of large fluctuations, the integral over earthquake ener- 
gies is dominated by the upper bound. The maximum 
energy E max (t) sampled by N(t)At earthquakes is given 

by the standard condition N(t)At J^ x{t) dE' P(E') k 
1. This yields the robust median estimate E max (t) ~ 
{N^At} 1 ^. Actually, E max (t) is itself distributed ac- 
cording to the Gutenberg-Richter distribution and thus 
exhibits large fluctuations from realization to realization, 
as we can see in Fig. 1. Putting this estimation of E max (t) 
in ©, we get 

N(t) cx £ dr {t _^l y+e IN(t)At}^/» . (4) 
Let us note the appearance of the new term 
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[N (r) At]*- 0- ^ resulting from the contribution of the 
upper bound in the integral J dE'P(E'). This term re- 
places the constant found for the case a < ji. Equation 
Q shows that the exploration of larger and larger events 
in the tail in the Gutenberg-Richer distribution trans- 
forms the linear Master equation into a non-linear 
equation: the non-linearity expresses a positive feedback 
according to which the larger is the rate N(t) of seismic- 
ity, the larger is the maximum sampled earthquake, and 
the larger is the number of daughters resulting from these 
extreme events. This process self- amplifies and leads to 
the announced finite-time singularity Q. However, to 
complete the derivation, we need to determine the yet 
unspecified time increment At. If N(t) obeys QJ, At 
is not a constant that can be factorized away: it is de- 
termined by the condition that, over At, N(t) does not 
change "significantly" in the interval [t, t + At], i.e., no 
more than by a constant factor. Using the assumed power 
law solution JJJ, this gives At oc t c — r. Using this and 
inserting Q into |J2J|, we get, 



" (a/ M )-l ' UC ^ 

(a/it) - 1 - e H(-9) 

m = ( n — i — > t c - 1 » c , 5 

[a/fx) - 1 

where H is the Heaviside function. Note that predicts 
an exponent to > 1 which is independent of 9 close to the 
critical time t c . This is due to the fact that the time de- 
cay of the Omori's kernel is not felt for t c — t < c, where c 
acts as an ultraviolet cut-off. It is also interesting to find 
that m = 1 independently of a and 9 in the regime 9 > 
(with of course a > fj) for which Omori's kernel ~ l/t 1+ 
decays sufficiently fast at long times that the predomi- 
nant contributions to the present seismic rate come from 
events in the immediate past of the present time of ob- 
servation. In constrast, the case 9 < is analogous to the 
anomalous long-time memory regime 

El 

which keeps for 
ever the impact of past events on future rates. 

This prediction, based on the careful analysis of the in- 
tegral in has been verified by direct numerical eval- 
uation of the equation We have also checked that 
numerical Monte Carlo simulations of the ETAS model 
generates catalogs of events following this prediction, in 
an ensemble or median sense. Figure I shows the cumu- 
lative number Af(t) — J Q dr N(t) of events for a typi- 
cal realization of the ETAS model and compares it with 
E m ax(t) to illustrate that Af(t) is mostly controlled by 
the sampling of E max (t), as discussed in the derivation of 
expression Q leading to the finite-time-singularity l[T]l. 
For the value li = 1 chosen here, E max {t) follows the same 
power law as the cumulative number, as observed. The 
dashed line is the power law prediction with JSJ for 
a/ li = 1.5 and 9 = —0.2 with slope to — 1 = 0.4. We have 
also generated 500 such catalogs and report in the inset 
the distribution d(m) of exponents to obtained by a best 



fit of M(t) for each of the 500 catalogs to a power law 
l/(ic — t) 7 ™ -1 - The median of d(m) is exactly equal to the 
prediction shown by the vertical thin line while the mode 
is very close to it. Note however a rather large dispersion 
which is expected from the highly intermittent dynamics 
characteristic of this extreme-dominated dynamics. We 
now report a few comparisons between the prediction (JSJ 
and the median value of the exponent to obtained from 
500 simulations for the following parameters: 9 = —0.2, 
a = 1.7, /i = 1, predicted to = 1.29, median to = 1.37; 
9 = —0.2, a = 1.3, /x = 1, predicted to = 1.67, median 
to = 1.61; 9 = —0.1, a = 1.5, ll = 1, predicted to = 1.20, 
median to = 1.29; 9 = —0.3, a = 1.5, fx = 1, predicted 
to = 1.60, median to = 1.62. For a > 1.8/z, the fluctu- 
ations are so large that a reliable determination of the 
median becomes questionable from a sample of 500 real- 
izations and many more would be needed. 

Figure 1 shows that the power law singularities are 
decorated by quite strong steps or oscillations, approxi- 
mately equidistant in the variable ln(i c — i). This log- 
periodicity has been previously proposed as a possi- 
bly important signature of rupture and earthquake se- 
quences approaching a critical point Q, Here, we 
present a simple novel mechanism for this observation, 
based on a refinement of the previous argument lead- 
ing to E max (t) - [iV(t)At]V^. Indeed, the most prob- 
able value for the energy E n of the n-th largest earth- 
quake ranked from the largest E\ = -E max to the small- 
est one is given by E n (t) = {[N(t)n + l]/[n/i + 1]} 1/m 
[3l|. where N{i) = j^N{t')dt'. Let us assume that 
the last new record was broken at time t\ leading to 
E x {ti) = {[Af(ti)fj, + l]/\p + 1]} 1/m . The next record 
will occur at a time > t\ whose typical value is such 
that E 2 (t2) = Ei(ti) (the last record Ei{t\) becomes the 
second largest event £2^2) when a new record Ei(t 2 ) oc- 
curs). For large Af(t), this gives %M = (2/i+l)/(/i+l). 
The prefered scaling ratio of the average log-periodicity 
is X = (t c -h)/(t c -t 2 ) - [(2 M +l)/(/i+ l)]V(m-i). For 
a/ li = 1.5,0 = —0.2, to = 1.4 corresponding to figure 1, 
we obtain A rj 2.3, which is compatible with the data. 

The prediction (JSJ) rationalizes the "inverse" Omori's 
law close to 1/ (t c — t ) t hat has been documented for 
earthquake forshocks [29j. The prediction © as well 
as the log-periodicity offers a general framework to ra- 
tionalize several previous experimental reports of precur- 
sory acoustic emissions rates prior to global failures [Ioj| . 
In this case, the energy release rate e(t) is found to fol- 
low a power law finite-time singularity. According to our 
theory, e(t) a N(t)E m ^(t) a l/(t c - £)-+(— D/f. 

Finally, we also show that this could explain star- 
quakes catalogs. Starquakes are assumed to be ruptures 
of a super-dense 1-km thick crust made of heavy nuclei 
stressed by super-strong stellar magnetic field. They are 
observed through the associated flashes of soft 7-rays ra- 
diated during the rupture. Starquakes exhibit all the 
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main stylized facts of their earthly siblings 30] . The 
thick line in figure 1 shows the cumulative number of star- 
quakes of the SGR1806-20 sequence, which is the longest 
sequence (of 111 events) ever attributed to the same neu- 
tron star, as a function of the logarithm of the time t c — t 
to failure. The starquake data is compatible with fj, = 1 
H3, a = 1.5 and 6 = -0.2, leading to m = 1.4. 

We are grateful to V. Keilis-Borok and V. Kossobokov 
for sharing the starquake data with us and W.-X. Zhou 
for discussions and help in a preliminary analysis of the 
data. 
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FIG. 1: Cumulative number of events (scale on the left) as 
a function of the time from the critical point t c for the star- 
quake sequence (solid black line) and one typical simulation 
of the ETAS model (solid thin line) generated with 9 — —0.2, 
a/n = 1.5 and c = 0.001. For the starquakes, t c is the time 
of the strongest observed starquake in the sequence. The 
dashed line shows the theoretical exponent m — 1 = 0.4 © 
for t a — t>c. The crosses x joined by straight segments give 
the time evolution of E max (t) (scale on the right). The inset 
gives the distribution of exponent measured for 500 numerical 
simulations. The median (vertical line) of the distribution of 
m-values is equal to the theoretical exponent m — 1.4 (for- 
mula JSJ). 
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